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Abstract 

Measurements of cosmic microwave background (CMB) anisotropy are ideal 
experiments for discovering the non-trivial global topology of the universe. To 
evaluate the CMB anisotropy in multiply-connected compact cosmological models, 
one needs to compute the eigenmodes of the Laplace-Beltrami operator. Using the 
direct boundary element method, we numerically obtain the low-lying eigenmodes 
on a compact hyperbolic 3-space called the Thurston manifold which is the second 
smallest in the known compact hyperbolic 3-manifolds. The computed eigenmodes 
are expanded in terms of eigenmodes on the unit three-dimensional pseudosphere. 
We numerically find that the expansion coefficients behave as Gaussian pseudo- 
random numbers for low-lying eigenmodes. The observed gaussianity in the CMB 
fluctuations can partially be attributed to the Gaussian pseudo-randomness of 
the expansion coefficients assuming that the Gaussian pseudo-randomness is the 
universal property of the compact hyperbolic spaces. 



1 Introduction 



In recent years, there has been a great interest in properties of CMB anisotropy 
in multiply-connected cosmological models |], 0, [3|, f|]. Most of these studies deal 
with flat models or non-compact hyperbolic models for which the eigenmodes are 
known explicitly Since no closed analytic expression of eigenmodes is known for 
compact hyperbolic (CH) models, so far, analysis of the CMB anisotropy in CH 
models has been considered to be quite difficult although they have interesting 
properties which are strikingly different from that of multiply-connected flat mod- 
els. For instance, in low Q Q adiabatic models, the large-angular fluctuations can be 
produced at periods after the last scattering as the curvature perturbations decay 
in the curvature dominant era. Therefore, the argument of the suppression of the 
large-angular fluctuations due to the "mode-cutoff" in the multiply-connected flat 
models cannot simply be applicable to the multiply-connected hyperbolic models. 

Because the effect of the multiply-connectedness becomes significant as the vol- 
ume of the space becomes small, it is very important to study whether the" small" 
universe scenario is plausible. For instance, the Weeks manifold and the Thurston 
manifold have volume~i? 3 where R denotes the curvature radius, and they are 
the smallest and the second smallest compact hyperbolic manifolds, respectively. 
For a technical reason, we study the properties of eigenmodes on the Thurston 
manifold. Since each type(scalar, vector and tensor) of perturbations evolves in- 
dependently on the locally homogeneous and isotropic FRW background space, in 
/c-space, the perturbations are given by the eigenmodes and the time evolution of 
the perturbations on the locally homogeneous and isotropic FRW space. The pe- 
riodic boundary conditions on the eigenmodes drastically change the nature of the 
CMB fluctuations on the topological identification scale while on smaller scale they 
asymptotically converge to that in the standard locally and globally homogeneous 
and isotropic FRW background space. Thus computation of eigenmodes are very 
important in understanding the properties of CMB fluctuations []. 

Computing the eigenmodes of the Laplace-Beltrami operator in CH spaces 

1 It should be noted that the computation of eigenmodes is also essential in the framework of 
spectral geometry ||. 
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(manifolds) is equivalent to solving the Helmholtz equation with appropriate peri- 
odic boundary conditions in the universal covering space. A number of numerical 
methods have been used for solving the Helmholtz equation such as the finite ele- 
ment methods and the finite difference methods || [7], || . 

A numerical method called the "direct boundary element method" (DBEM) has 
been used by Aurich and Steiner for finding out the eigenmodes of the Laplace- 
Beltrami operator in a two-dimensional compact multiply-connected space for 
studying the statistical properties of the eigenmodes in highly-excited states or 
equivalently the semi-classical wavefunctions || . We find that this pioneering work 
for "quantum chaology" , the study of the imprints of classical chaos in the quan- 
tum mechanical counterparts is very useful for the study of the CMB anisotropy 
in CH cosmological models as well. The advantage of the DBEM is that it reduces 
the dimensionality of the problem by one which leads to economy in the numer- 
ical task. Since one needs to discretize only the boundary, generation of meshes 
is much easier than the other methods. Furthermore, as we shall see later, the 
DBEM is suitable for expanding the eigenfunctions that are continued onto the 
whole Poincare ball by the periodic boundary conditions in terms of eigenfunc- 
tions on the pseudosphereQ. In order to compute a set of expansion coefficients, 
one needs to compute the values of the corresponding eigenfunction on a hyper- 
bolic sphere with appropriate radius. If one does not care about the normalization 
of the eigenfunctions, unlike the FEM, the computation of the eigenfunctions on 
the whole fundamental domain is not necessary. In the DBEM, computation of 
the eigenfunction at an arbitrary point either inside or outside of the fundamen- 
tal domain can be done by using the values of the eigenfunction and the normal 
derivatives on the boundary. 

As the classical dynamical systems in CH spaces are strongly chaotic, one can 
naturally assume that the imprint of the classical chaos is hidden in the corre- 
sponding quantum systems in some way. It has been found that the expansion 
coefficients with a certain basis behave as if they are random Gaussian numbers 



in some classically chaotic systems TU|), which is consistent with the predic- 



tion of random-matrix theory([TT|). Since the CMB temperature fluctuations in 

2 A set of the continued eigenfunctions is a subset of all eigenfunctions on the universal covering 
space. 
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CH spaces are written in terms of expansion coefficients and the eigenfunctions 
on the universal covering space plus initial fluctuations, if the random behavior of 
the expansion coefficients is confirmed, the origin of the random gaussianity in the 
CMB temperature fluctuations can be partially explained in terms of the geometric 
property of the universe. 

In this paper, we introduce the DBEM for solving the Helmholtz equation. 
Then we apply the DBEM for computing the low-lying eigenmodes of the Laplace- 
Beltrami operator in the Thurston manifold. The computed eigenfunctions are 
naturally continued onto the whole Poincare ball because of the periodic boundary 
conditions and are expanded in terms of eigenmodes on the simply-connected pseu- 
dosphere. Statistical properties of the expansion coefficients are examined, since 
they are key factors in understanding of the CMB anisotropy in CH models. 

2 The direct boundary element method (DBEM) 

The boundary element methods (BEM) use free Green's function as the weighted 
function, and the Helmholtz equation is rewritten as an integral equation defined 
on the boundary using Green's theorem. Discretization of the boundary integral 
equation yields a system of linear equations. Since one needs the discretiztion on 
only the boundary, BEM reduces the dimensionality of the problem by one which 
leads to economy in the numerical task. To locate an eig envalue, the DBEM Q re- 
quires one to compute many determinants of the corresponding boundary matrices 
which are dependent on the wavenumber k. 

Firstly, let us consider the Helmholtz equation with certain boundary condi- 
tions, 

(V 2 + fc 2 )u(x) = 0, (1) 

which is defined on a bounded M-dimensional connected and simply-connected 

domain Q which is a subspace of a M-dimensional Riemannian manifold M. and 

the boundary dfl is piecewise smooth. V 2 = V l Vj, (i = 1,2,-- -, M), and Vj is the 

3 The DBEM uses only boundary points in evaluating the integrand in Eq.(||). The indirect 
methods use internal points in evaluating the integrand in Eq.rtq) as well as the boundary points. 
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covariant derivative operator defined on M.. A function u in Sobolev space H 2 (Q) 
is the solution of the Helmholtz equation if and only if 

K[u(x),v(x)] = ((V 2 + k 2 ) tz(x),«(x)\ = 0, 



(2) 

where t> is an arbitrary function in Sobolev space H 1 (fl) called weighted function 
and ( ) is defined as 

(a, 6) = / ab^dV. (3) 
Next, we put u(x) into the form 



M 
3=1 



(4) 



where 0j's are linearly independent square-integrable functions. Numerical meth- 
ods such as the finite element methods try to minimize the residue function 1Z for 
a fixed weighted function u(x) by changing the coefficients Uj. In these methods, 
one must resort to the variational principle to find the u/s which minimize 1Z. 

Now we formulate the DBEM which is a version of BEMs. Here we search 
w(x)'s for the space C 1 ^) fl C 2 (Q) R L 2 (fi). First, we slightly modify Eq.© using 
the Green's theorem 



(V 2 u)vy/gdV - / (V^w^c^/ 



{Viv)u^gdS\ (5) 



'n Jn Jan Jan 

where g = det {#,,,•} and rfV" = rfxi . . . dxu\ the surface element dS l is given by 

1 
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(6) 



where e 



31—3M+1 



denotes the M+l-dimensional Levi-Civita tensor. Then Eq.(EI) be- 



comes 



(V 2 v + k 2 v)uy/g dV + / (Wiu)vy/gdS l - / (Viv)u^gdS 1 = 0. (7) 
<> Jan Jan 



As the weighted function v, we choose the fundamental solution Ge(x, y) which 
satisfies 

(V 2 + E)G E (x,y) = 6 D (x-y), (8) 

where E = k 2 , and 5d(x — y) is Dirac's delta function. G E (x, y) is also known as 
the free Green's function whose boundary condition is given as 

lim G E (x,y) = 0, (9) 

<2(x,y)-^oo 

where d(x, y) is the geodesic distance between x and y. Let y be an internal point 
of Q. Then we obtain from Eq.(^) and Eq.(||), 

u(y)+ [ G E (x,y)V&JgdS i - [ (ViG E (x,y))u^g dS' = 0. (10) 
Jan Jan 

Thus the values of eigenfunctions at internal points can be computed using only 

the boundary integral. If y G dQ, we have to evaluate the limit of the boundary 

integral terms as Ge(x, y) becomes divergent at x = y (see appendix A). The 

boundary integral equation is finally written as 

\<y)+ I GEix^ViUy/gdS 1 - f (V.G^x, y))u ^dS l = 0, (11) 
1 Jan Jan 

or in another form, 



\u{y)+l G E (x,y)^n l ^dS- [ d ° E ^ y) n'ti ^ dS = 0, (12) 
1 Jan ox Jan ° x 



where n l = dS l /dS and dS = \f~dl¥d~Si. Note that we assumed that the boundary 
surface at y is sufficiently smooth. If the boundary is not smooth, one must 
calculate the internal solid angle at y (see appendix A). Another approach is to 
rewrite Eq.(|T0|) in a regularized form |fT5| . We see from Eq. (|TTD or Eq. (|T2"D that 
the approximated solutions can be obtained without resorting to the variational 



principle. Since it is virtually impossible to solve Eq.(12) analytically, we discretize 
it using boundary elements. Let the number of the elements be N. We approximate 
u by some low-order polynomials (shape function) on each element as u = C\ + 
C2 V + C3 £ where 77 and £ denote the coordinates on the corresponding standard 
element 0. 

4 It can be proved that the approximated polynomial solutions converge to u(x) as the number 
of boundary elements becomes large jl4|, |l5) . 
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Then we have the following equation: 

[H]{u} = [G]{q}, q=?j±, (13) 

where {u} and {q} are N-dimensional vectors which consist of the boundary val- 
ues of an eigenfunction and its normal derivatives, respectively. [H] and [G] are 
N x N- dimensional coefficient matrices which are obtained from integration of the 
fundamental solution G\e(x, y) and its normal derivatives multiplied by and 
respectively. The explicit form of [H] and [G] for constant elements are given in sec- 
tion 4. Note that the elements in [H] and [G] include k implicitly. Because Eq.flTB"!) 
includes both u and q, the boundary element method can naturally incorporate 
the periodic boundary conditions: 

u(x) =u(0i(x)), g(x) = -g(^(x)), on dSl, (14) 

where g^s are the face-to-face identification maps defined on the boundary (see 
appendix B). The boundary conditions constrain the number of unknown constants 
to N. Application of the boundary condition ( |T4"D to Eg. ([IB]) and permutation of 
the columns of the components yields 

[A}{x} = 0, (15) 

where NxN- dimensional matrix A is constructed from Gij and and N- dimensional 
vector x is constructed from Wj's and q^s. For the presence of the non-trivial solu- 
tion, the following relation must hold, 

detL4] = 0. (16) 

Thus the eigenvalues of the Laplace-Beltrami operator acting on the space C 1 (Cl) n 
C 2 (Q) R L 2 (Q) are obtained by searching for fe's which satisfy Eq . (|16"D . 
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3 Computation of low-lying eigenmodes 



In this section, we apply the DBEM for computing the low-lying eigenmodes 
on the Thurston manifold Q2- We have chosen Q2 for a technical reason that 
the fundamental domain of Q2 is much simpler than that of the Weeks manifold 
Qi (generation of meshes is much simpler). See appendix B for understanding the 
basic aspects of three-dimensional hyperbolic geometry. 

The Helmholtz equation in the Poincare coordinates is written as 



1 



x - V, 



u + k 2 u = 0, 



;i7) 



where and Ve are the Laplacian and the gradient on the corresponding three- 
dimensional Euclidean space, respectively. Note that we have set the curvature 
radius R = 1 without loss of generality. By using the DBEM, the Helmholtz 
equation (|17]) is converted to an integral representation on the boundary. Here 



Eq.([12]) can be written in terms of Euclidean quantities as 



My) + 



where dS = 2(1 



G k {x,y)—n E dS 

an ox 1 



<9G fc (x,y) 



un l rp dS = 0, 



an dx i 

2 ) _1 dS^. The fundamental solution is given as WB, [T7 



Gfc(x ' y) = ~47 



a 



+ vV 5 



71 71 

- < arg s < -, 



(19) 



where s = yl — k 2 and a = cosh d(x, y) 
boundary elements T j as 

N r 



Then Eq.([l8|) is discretized on the 



1 



E 

j=i 



Gfc(x/,yj) dn dS 



dG k (x!,yj] 
dn 



u(yj) dS 



0, (20) 



where N denotes the number of the boundary elements. An example of iV = 1168 
elements on the boundary of the fundamental domain in the Poincare coordinates 
is shown in figure [l]. These elements are firstly generated in Klein coordinates in 
which the mesh-generation is convenient. The maximum length of the edge Al in 
these elements is 0.14. The condition that the corresponding de Broglie wavelength 
27r/fc is longer than the four times of the interval of the boundary elements yields 
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Figure 1: 1168 boundary elements 
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a rough estimate of the validity condition of the calculation as k<ll. On each Tj, 
u and q = du/dn are approximated by low order polynomials. For simplicity, we 
use constant elements: 



u 



xj) = u J = Const. , <z(xj) = q J = Const. , on Tj. (21) 



Substituting Eq.(pTl) into Eq. (p0|) , we obtain 

TV N 



j=i j=i 



HlJ = \ TT I 



iJ/ j — 2 I = J, 



H u = I ^( Xl , yj )dS(yj), G u = I G k (xi,yj)dS(yj). C22) 



where 



The singular integration must be carried out for I-I components as the fundamental 
solution diverges at (xj = yi). This is not an intractable problem. Several nu- 



merical techniques have already been proposed by some authors |T3|, |TU| . We have 
applied Hayami's method to the evaluation of the singular integrals [19|]. Introduc- 
ing coordinates similar to spherical coordinates centered at x/, the singularity is 
canceled out by the Jacobian which makes the integral regular. 

Let gi [i = 1, 2, . . . , 8) be the generators of the discrete group T which identify a 
boundary face Fi with another boundary face gi(Fi): 

gi(xi) = Xj, x.j e Fi. (23) 

The boundary of the fundamental domain can be divided into two regions 8Qa 
and 8Qb and each of them consists of N/2 boundary elements, 

dn A = UF u dn B = Ug i (F i ), t = l,2,..., 8. (24) 

The periodic boundary conditions 

u(gi(xi)) = ufa), q(gi(xi)) = -q(xi), i = l,2,...,8 (25) 
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reduce the number of the independent variables to N, i.e. for all xg G OQb, there 
exist Qi G T and G dfi^ such that 

u(x B ) = «(^(x A )) = u(x A ), g(x fl ) = -g(flfi(xx)) = -<?(xa). (26) 

Substituting the above relation into Eq. ([22]) , we obtain 

Haa + Hab —Gaa + Cms 
-Hba + Hbb —Gba + 

where Ua = (u 1 , u 2 , . . . u N l 2 ) and qA = (q 1 , q 2 , ■ ■ ■ q N ^ 2 ) and matrices H = {Hjj} 
and G = {Gu} are written as 

jj _ Haa H A b 
[ H B a H B b 

Eq. (p7[) takes the form 

[A(A;)]{:r} = 0, (29) 

where Nx N- dimensional matrix A is constructed from G and H and N- dimensional 
vector x is constructed from and qA- For the presence of the non-trivial solution, 
the following relation must hold, 

det{A(k)] = 0. (30) 

Thus the eigenvalues of the Laplace-Beltrami operator in a CH space are obtained 
by searching for fc's which satisfy Eq.(^U]). In practice, Eq.([3(]) cannot be exactly 
satisfied as the test function which has a locally polynomial behavior is slightly 
deviated from the exact eigenfunction. Instead, one must search for the local 
minima of det[A(k)]. This process needs long computation time as A(k) depends 
on k implicitly. Our numerical result (k< 10) is shown in table |l[ 

The first "excited state" which corresponds to k = k\ is important for the 
understanding of CMB anisotropy. Our numerical result ki = 5.41 is consistent 
with the value 5.04 obtained from Weyl's asymptotic formula 

m= y°m* t „ s ^3i, „»i, 

07T Z 

assuming that no degeneracy occurs. One can interpret the first excited state 
as the mode that has the maximum de Broglie wavelength 1njk\. Because of 



Qa 



(27) 



G 



Gaa Gab 
Gba Gbb 



(2i 
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k 


m k 


5.41 


1 


5.79 


1 


6.81 


1 


6.89 


1 


7.12 


1 


7.69 


1 


8.30 


1 


8.60 


1 


8.73 


1 


9.26 


2 


9.76 


1 


9.91 


1 


9.99 


1 



Table 1: Eigenvalue k and multiplicity 



the periodic boundary conditions, the de Broglie wavelength can be approximated 
by the "average diameter" of the fundamental domain defined as a sum of the 
inradius r_ and the outradius which yields k\ = 4.9 just 10% less than the 
numerical value. From these estimates, supercurvature modes in small CH spaces 
(Vol(A4) ~ 1) are unlikely to be observed. 

To compute the value of eigenfunctions inside the fundamental domain, one 
needs to solve Eq.(^). The singular decomposition method is the most suitable 
numerical method for solving any linear equation with a singular matrix A, which 
can be decomposed as 

A = U^DV, (32) 

where U and V are unitary matrices and D is a diagonal matrix. If Da in D is 
almost zero then the complex conjugate of the i-th row in V is an approximated so- 
lution of Eq.(^). The number of the "almost zero" diagonal elements in D is equal 
to the multiplicity number. Substituting the values of the eigenfunctions and their 
normal derivatives on the boundary into Eq. (fl~0|) , the values of the eigenfunctions 
inside the fundamental domain can be computed. Eigenfunctions k = 5.41 and 

5 The inradius r_ is the radius of the largest simply-connected sphere in the fundamental 
domain, and the outradius r + is the radius of the smallest sphere that can enclose the fundamental 
domain. ?*_ =0.535, r + = 0.7485 for the Thurston manifold. 
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k=5 .41 



k= 9 . 9 9 




A; = 9.99 in Poincare coordinates plotted as (xi,X2,h), where h = u v (xi, X2, 0.0087) 
are shown in figure 0. The eigenfuctions we computed are all real-valued. Note 
that the non-degenerated eigenfunctions must be real-valued. 

The numerical accuracy of the obtained eigenvalues is roughly estimated as 
follows. First, let us write the obtained numerical solution in terms of the exact 
solution as k = ko + 8k and Wfc(x) = u ko (x) + 8u k (x), where k and u ko (x) are 
the exact eigenvalue and eigenfunction, respectively. The singular decomposition 
method enables us to find the best approximated solution which satisfies 

[A]{x} = e, |e|«l, (33) 

where e is a N-dimensional vector and | | denotes the Euclidean norm. It is expected 
that the better approximation gives the smaller |e|. Then Eq. (1331) can be written 



/ G fco+a (x,yj)(A + (ko + 8kf)(u ko (x) + 5u k (x)) ^gdV x = e( yj ). (34) 
Jn 

Ignoring the terms in second order, Eq. (|3^) is reduced to 

/ G fc (x, yj )((A + kl)8u k (-K) + 2k8ku k {x)) ^gdV x = e(yj). (35) 
Jn 

Since it is not unlikely that (A + kQ)8u k (x) is anticorrelated to 2k8ku k {x), we 
obtain the following relation by averaging over y j, 

2k\6k\(\J^G k {x,yj)u k {x)y/gdvJ^ ~ < |e| >, (36) 

13 



1 

0.5 



. 1 
. 05 



.01 
D .005 



V 6 .0001 . 01 



1 

0.5 



0.1 
. 05 



.01 
.005 



*•. Tr «*«" 



n -6 0.0001 0.01 



Figure 3: Error (err) versus hyperbolic distance to the boundary for 291 points in- 
side the fundaental domain, d denotes the hyperbolic distance from the evaluation 
point x to the nearest point on the boundary. 



where <> denotes the averaging over yj Thus one can estimate the expected 
deviation of the calculated eigenvalue \Sk\ from Wfc(x) and e(yj). We numerically 
find that \6k\ = 0.005 for k = 5.41 and \Sk\ = 0.01 for k = 9.91. The other 
deviation values lie in between 0.005 and 0.01. 

By computing the second derivatives, one can also estimate the accuracy of the 
computed eigenfunctions. The accuracy parameter err is defined as 

err(£;,x) = (A + k 2 )u k (x), (37) 

where ttfc(x) is normalized (0(«fc(x)) ~ 1). We see from figure |3] that the accu- 
racy becomes worse as the evaluation point approaches the boundary. However, 
for points with hyperbolic distance d > 0.1 between the evaluating point and the 
nearest boundary, the errors are very small indeed: err < 10~ 4 ~~ 5 . This result is 
considered to be natural because the characteristic scale L of the boundary ele- 
ments is ~ 0.07 for our 1168 elements. If d < L, the integrands in Eq.fllOD become 
appreciable on the neighborhood of the nearest boundary point because the free 
Green's function approximately diverges on the point. In this case, the effect of 
the deviation from the exact eigenfunction is significant. If d » L, the integrand 
on all the boundary points contributes almost equally to the integration so that 
the local deviations are cancelled out. 

As we shall see in the next section, expansion coefficients are calculated using 
the values of eigenfunctions on a sphere. Since the number of evaluating points 
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which are very close to the boundary is negligible on the sphere, expansion coeffi- 
cients can be computed with relatively high accuracy. 
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4 Statistical properties of eigenmodes 



Properties of eigenmodes of the Laplace-Beltrami operator are determined by 
the Helmholtz equation. Therefore, at first glance it does not seem to make a sense 
to study the statistical properties of the eigenmodes. 

However, if one recognizes the Laplace-Beltrami operator in a CH space as the 
Hamiltonian in a quantum mechanical system, each eigenmode can be interpreted 
as a wavefunction in an eigenstate. Since the corresponding classical system is 
known to be a typical chaotic system (K-system), it is natural to assume that the 
imprints of classical chaos is hidden in the corresponding quantum system. 

Recent studies have demonstrated that some of the statistical properties of 
energy spectrum are in accordance with the universal prediction of random-matrix 



theory(RMT) pO], |2T|. RMT also predicts that the squared expansion coefficients 
| Qi | 2 of an eigenstate with respect to a generic basis are Gaussian distributed R [TO, 



TTJ. In the limit iV — > oo, x = |aj| 2 obeys the statistics given for three universality 
classes of the orthogonal (GOE, /i = 1), unitary (GUE, /z = 2) and symplectic 
(GSE, fi = 4) ensembles, each distribution function P is given by 

P {x) = (^X' rO/2) x^' 2 - 1 e~^ 2 . (38) 



% 2, 

In our case, as the time-reversal symmetry of the Hamiltonian implies, one expects 
that \a,i\ 2 obeys the GOE prediction. In order to apply the GOE prediction to 
the statistical properties of eigenstates on CH spaces, one needs to find a set of 
orthonormal basis but no closed analytic expression is known for any CH spaces. 
To avoid the problem, Aurich and Steiner noticed that the wavefunctions on the 
hyperbolic octagons can be continued onto the universal covering space H 2 , and 
eigenstates can be expanded in terms of circular-waves [[J. They numerically found 
that the squared expansion coefficients obeys the GOE prediction in highly excited 
states of a hyperbolic asymmetrical octagon model. 

We extend their method to three-dimensional CH models where we consider 
only low- lying modes. First, we normalize the obtained 14 eigenf unctions on the 
Thurston manifold. The eigenfunctions are naturally continued onto the whole 
unit Poincare ball by the periodic boundary condition. As a "generic basis", we 
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consider a set of orthonormal eigenfunctions Q v i m (T-complete functions) on the 
unit pseudosphere which is isometric to the Poincare ball, 

Qulm = X„i(x)Y lm (9,(f)), 



X vl ee (-l)' +1 ^(ll(^ + ^ 2 )) ^sinh^ 



d(coshx) l+1 



r(nH«) ./Xp-»/ VA ✓ = (39) 



r(i/i) V sinhx 1/4-1/2 
where P, Yj m and V denote the associated Legendre function, the spherical har- 
monics and gamma function, respectively. P can be written in terms of the hyper- 



geometric function 2 JF X [22 



p-i-y 2 2 (cosh X ) = (C °^! X) t 1/2 *Fi i\ - *> £ + + s-h 2 i X ) . (40) 



ui-i/2 y—Aj - r( | + z) ^ 2 ' ' 2 ' ' 2 

Eigenfunctions ttj, can be expanded in terms of Q vim's as 

iu = ^^imX vl { X )Y lm {9,(j)). (41) 



Ml 

lm 



Note that each has no components with u ^ v because Q^ m 's are complete 
and linearly independent. 

At first glance the computation of in Eq. ([41]) seems cumbersome as the 
domain of the integration extends over the whole pseudosphere. Fortunately, one 
can obtain ^im by evaluating two-dimensional integrals. £ u i m can be written as 

UmX vl {xo) = [ u„( Xo , o, 0) Y lm (e, <f>)dn, (42) 



which is satisfied for the arbitrary value of Xo- in practice the numerical instability 
occurs in the region where the absolute value of X„i is too small. In our computa- 
tion, the values of Xo are chosen as shown in table ||. Thus can be computed 
if one obtains the values of eigenfunctions on the sphere 

x\ = tanh — sin 6 cos <p, x 2 = tanh — sin 6 sin <p, x 3 = tanh — cos 6 (43) 
2 2 2 

with radius Xo- 

In order to compute the values of eigenfunctions on the sphere with radius longer 
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<1<8 


8 < 1 < 13 


13 < 1 < 18 


k<8 


0.53 


1.3 


1.6 


k>8 


0.53 


1.1 


1.3 



Table 2: An example of choice of Xo for which the absolute value of X v i{xo) is not 
too small. 



than the inradius r_ =0.535, the points outside of the fundamental domain must 
be pulled back to the inside, since Eq.(|H]) is valid only if y is a set of coordinates 
of an internal point. 

The plots of eigenfunctions on a sphere Xo = 1-6 are shown in figure (| and figure 
|5|. Apparent structure of the eigenfunctions on the sphere seems complicated. 
However, some regular patterns are hidden in the structure due to the periodic 
boundary conditions. Actually, there are pairs of highly correlated points on the 
sphere, since any partial surface Su of the sphere that is enclosed by copies of 
the boundary of the fundamental domain pulled back to inside the fundamental 
domain by the corresponding element of the discrete isometry group intersects 
another partial surface S 1 ^ that is pulled back to inside the fundamental domain. 
To evaluate the correlation pattern, let us estimate how often a sphere with radius 
X = Xo intersects the copies of the fundamental domain. The approximate number 
rii of the copies of the fundamental domain inside the sphere with radius (in proper 
length) Xo is given by 

_ 7r(sinh(2x ) - 2 Xo ) . . 

Vol(Q 2 ) ■ l44j 
From this formula, in the case of the Thurston manifold, ni ~ 29 if Xo = 1-6. 

Because the sphere intersects the fundamental domain at random, the copies of 

the fundamental domain on the sphere stick out their half portions on average. 

Therefore, the approximate number n 2 of the copies that intersect the sphere is 

given by 

_ 7r(sinh(2(x + r ave )) - sinh(2(x D - r ave )) - Ar ave ) 

112 ~ v^m ' (45) 

where r ave = (r + +r_)/2. This estimate gives n 2 ~120 if x G = 1.6. Approximating 
each eigenmode by de Broglie waves, we obtain the corresponding fluctuation scale 
SA in steradian on the sphere, 

47T 47T 2 

SA = — 

ri2 k 



16vr 2 Vol(Q 2 ) , \ 

(46) 



fc 2 (sinh(2(x + r ave )) - sinh(2(x o Tave 

)) -4r a „ e )' 

Thus correlation patterns are observed in pairs of patches with typical size 8A. 
When Xo = 1-6, angular fluctuation scales are given as 51 ~ \^5A = 21°, 12°, for 
k = 5.41, 9.99, respectively. 

Next, we extract a set of independent variables from £„/ m 's. In general, any 
Q u i m is related to Q u i- m as 

Qvi-m=(-l) m Qli m F(v,l), (47) 

where 

If u v is real, from Eq. (f4"7j) , 

^ ^ ^ulm Qulm 
I m 



= \^ f* O* 

/ j ^ulm vim 

I m 

= ^(-^Cl-rnQvlmF-^l), (49) 
I rn 

therefore, 

U m = {-l) m Ci- m F-\^l). (50) 

Thus ^z_ m can be written in terms of Culm- To extract a set of independent 
variables from ^ m 's, we rewrite Eq. (f49|) as follows 



U v ^ ^ ^!/im Qulm 4" ^ ^ £i4m Qulm 

l,m<0 l,m>0 

= -Im(^oo) Im^oo) + V^i^^iojR. 



ulO 

l>0 



l>0,m>0 

where 



+ v / 2(Re(e,/ m )v / 2Re(g^ m ) - lm{U m )V2lm{Q ulm )j 



(51) 



R,io = {cui)- 112 ( Re(Qwo) + Yrf&iJ Im(g ^ o)i ' ' (52) 
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Figure 4: Absolute values of eigenfunction u(k = 5.41) on sphere Xo = 1-6 
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Figure 5: Absolute values of eigenfunction u(k = 9.99) on sphere Xo = 1-6 
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and 



C U 1 



x Im 



2F(u,l) 

T(l + vi + l) 
T(vi) 



1-F(v,l) \* 
1 + F(v,l) 



+ 1 



Vim) 



T(l + vi + l) 



(53) 



l + Re(F(u,l)Y 

Thus the real eigenfuctions can be written in terms of real independent coefficients 
a u [ m and real- valued R u i m , 

U v = a u lm Rylmi (54) 



where 



l, m 



and 



«i/oo = — im (^oo), = v / ^R- e (^o), I > 0, 

a„j m = V2R&(£, v i m ), m > 0, a„ Zm = -\/2Im(^_ m ), m < 0, (55) 

-Ri/00 = im (Q^Oo); 

R v im = \p2Re(Q u i m ), m > 0, i?^ m = v2Im(Q^_ m ), m < 0. (56) 



Now we turn to the statistical properties of the coefficients a v \ m . As in 
consider the cumulative distribution of following quantities, 

I O'ulra Q>v I 



(57) 



where a u is the mean of a^ m 's and a 2 v is the variance. The cumulative distribution is 
compared to the cumulative RMT distribution functions which are directly derived 
from Eq.pl), 



I^(x) = / dx'P^x) 
Jo 

= 7( ^ 2 'f/ 2) , (58) 

where 7(2;, y) is the incomplete gamma function. To test the goodness of fit between 
the computed cumulative distribution function and that predicted by RMT, we use 
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Kolmogorov-Smirnov statistic which is the least upper bound of all pointwise 
differences \In(x) — I(x)\ f23f , 



Djv = sup \In(x) — I(x)\, (59) 

X 

where In{x) is the empirical cumulative distribution function defined by 

{0, x < 

j/N, y j <x<y j+1 , j = 1, 2, . . . , N-l, (60) 
1, Vj < x, 

where y\ < y 2 < ■ ■ ■ < Vn are the computed values of a random sample which 
consists of N elements. If In(x) is "close" to I(x), the observed D N must be so 
small that it falls within the range of possible fluctuations which depend on the 
size of the random sample. For the random variable for any z > 0, it can be 
shown that the probability of < d is given by [fj4[] 

lim P(D N < d = zN- 1/2 ) =L(z), (61) 

where 

oo 

L{z) = 1-2 ^(-ly-V 2 -^ 2 . (62) 

3=1 

From observed maximum difference = d, we obtain the significant level ao = 
1 — P which is equal to the probability of D N > d. If od is found to be large 
enough, the hypothesis In(x) = I(x) is verified. The computed cumulative dis- 
tributions of |a^ m | 2 and the GOE(/x = 1) prediction h(x) for four examples are 
plotted in figure || and the maximum difference d and the significant levels a d for 
< / < 10 and < / < 18 are shown in table |3]. Note that the last digit in is 
not guaranteed, since Eq.(pH) is an asymptotic formula. 

We see from figure |6] and table [5] that the agreement with GOE prediction is 
remarkably good. The Gaussian behavior for highly excited states is naturally 
expected as the semiclassical wavefunctions must reflect the chaotic nature in the 
corresponding classical systems. However, the Gaussian behavior for low-lying 
modes is not apparent at all. It is possible that the non-gaussianity behavior is 
rather prominent as these modes have fluctuations on large scale and that reflects 
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Figure 6: The cumulative distribution of \a vim — a u i m \ 2 /a 2 are plotted with the 
theoretical curve for GOE prediction (full curve). The statistics are based on 121 
expansion coefficients for I < 11 and 361 coefficients for I < 19. 



the pure quantum mechanical behavior^. Nevertheless, our numerical results serve 

to strengthen the hypothesis that the expansion coefficients behave as Gaussian 

pseudo-random numbers even for low-lying modes. 

6 The typical scale of angular fluctuation of I mode is approximately it/ (21) and the typical 
scale of radial fluctuation of k mode is approximately given by 2ir/k. 
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< I <10 


< I <18 


k 


D N x 10 2 


OLD 


D N x 10 2 




5.41 


9.42 


23.4% 


4.29 


51.8% 


5.79 


3.88 


99.3% 


4.29 


52.1% 


6.81 


5.95 


78.6% 


3.47 


77.8% 


6.89 


9.34 


24.1% 


2.36 


98.8% 


7.12 


7.12 


57.2% 


2.59 


96.8% 


7.69 


10.23 


15.9% 


4.84 


36.7% 


8.30 


6.38 


70.8% 


3.88 


65.0% 


8.60 


5.63 


83.8% 


2.09 


99.7% 


8.73 


9.46 


22.9% 


2.78 


94.3% 


9.26a 


7.21 


55.6% 


3.46 


78.1% 


9.26b 


5.99 


77.9% 


6.11 


13.5% 


9.76 


7.41 


52.0% 


4.57 


43.8% 


9.91 


8.90 


29.3% 


3.23 


84.7% 


9.99 


4.27 


98.0% 


3.72 


70.0% 


ave. 


7.23 


56.3% 


3.69 


68.8% 



Table 3: The Kolmogolov-Smirnov statistics D N and the significance levels ao for 
the test of the hypothesis In(x) ^ I(x) and their averages. N = 121 for < / < 10 
and iV = 361 for < I < 18. The mode k = 9.26 is degenerated into two modes, 
which (after orthogonalization)are denoted by k = 9.26a and k = 9.266. 
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Next we examine the randomness of a^ m 's. Because a„z m 's are actually deter- 
mined by the Helmholtz equation, it is appropriate to describe a^ m 's as pseudo- 
random numbers. We apply the run test for testing randomness (see |23|). 



Suppose that we have n observations of the random variable X and m obser- 
vations of the random variable Y. The combination of those variables into m + n 
observations placed in ascending order of magnitude yields 

xxx yy xx yyy x y xx yy, 

where x denotes an observation of X and y denotes an observation of Y. Each 
underlined group which consists of successive values of X or Y is called run. The 
statistics of number of runs are used for testing whether X and Y have the same the 
distribution function. Regardless of the type of the distribution function, the run 
number r is known to behave as Gaussian random numbers in the limit m, n — > oo. 

The run test is also used as a test for randomness. Let a\, 0,2, ■ ■ ■ , ajy be the 
observed values of a random variable A. For simplicity, assume that N is even. The 
median divides the observed values into a lower and an upper half. It is represented 
as L if it falls below the median, and it is represented as U if it falls above the 
median. For instance, a sequence 

UUU L UU LLUUL UU, 

has 8 numbers of runs (r = 8). The critical region for testing the hypothesis of 
randomness is of the form r < C\ or r > C2 where C\ and C2 is readily given by the 
Gaussian distribution function. The significant level a r is the probability of r > c\ 
or r > C2. As the Kolmogolov-Smirnov test, a r is given by the observed r. 

The run numbers r and the significant levels a r are shown in table f|. High 
significant levels are again observed except for the one at /c = 8.73 for < / < 18. 
As the corresponding r is larger than the averaged value, this may be due to the 
cyclic effect. On the whole, it is concluded that a v \ m 's behave as if they are random 
variables. 
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< / <10 


< I <18 


k 


r 


a r 


r 


a r 


5.41 


62 


85.5% 


185 


67.3% 


5.79 


58 


46.5% 


174 


39.9% 


6.81 


69 


14.4% 


196 


11.4% 


6.89 


60 


71.5% 


168 


14.0% 


7.12 


69 


14.4% 


184 


75.2% 


7.69 


57 


36.1% 


191 


29.2% 


8.30 


63 


71.5% 


177 


59.8% 


8.60 


59 


58.3% 


184 


75.2% 


8.73 


70 


10.0% 


201 


3.5% 


9.26a 


56 


27.3% 


177 


59.8% 


9.26b 


55 


20.1% 


182 


91.6% 


9.76 


59 


58.3% 


182 


91.6% 


9.91 


70 


10.0% 


196 


11.4% 


9.99 


58 


46.5% 


179 


75.2% 


ave. 


61.8 


40.7% 


184 


50.4% 



Table 4: The run numbers r and the significance levels a r for the test of the 
hypothesis that a u i m 's are random variables. N = 121 for < I < 10 and N = 361 
for < I < 18. The mode k = 9.26 is degenerated into two modes, which (after 
orthogonalization)are denoted by k = 9.26a and A; = 9.266. 
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5 Summary 



In this paper, we have demonstrated that the DBEM is eminently suitable for 
computing eigenmodes on CH spaces and we obtain some low-lying eigenmodes on 
a CH space called Thurston manifold which is the second smallest in the known 
CH manifolds and we have studied the statistical properties of these eigenmodes. 

The low-lying eigenmodes are expanded in terms of eigenmodes on the pseudo- 
sphere, and we find that the expansion coefficients behave as if they are Gaussian 
random numbers. Why are they so random even for low-lying modes? It should 
be pointed out that the randomness of the expansion coefficients for low-lying 
eigenmodes is not the property of the eigenmodes themselves but rather the prop- 
erty of the images of eigenmodes on the whole universal covering space, since the 
fluctuation scales for low-lying eigenmodes are comparable to the the size of the 
fundamental domain. We conjecture that the origin of the random behavior of 
eigenmodes comes from the almost randomly distributed images of a set of points 
in the universal covering space. 

Computation of eigenmodes is essential in simulating the CMB in CH cosmo- 
logical models. As the DBEM needs only a set of face-to-face identification maps 
and the discretization of the corresponding fundamental domain, it can be applied 
to other CH spaces straightforwardly. However, the computation of the modes 
with small fluctuation scale k » 1 is still a difficult task as the number of modes 
increases as N cx k 3 . 

Nevertheless, the contribution of the modes with small fluctuation to the tem- 
perature correlation of CMB can be estimated by assuming that the expansion 
coefficients for excited states (k » 1) also behave as Gaussian pseudo-random 
numbers as well as that for low-lying modes. The assumption is numerically con- 
firmed in a two-dimensional CH model |J. 

If the observed Gaussian pseudo-randomness is found to be the universal be- 
havior in CH spaces for low-lying modes as well as excited modes, the origin of the 
gaussianity of the CMB fluctuations can be partially explained. This is because 
the amplitude of the CMB fluctuation is written in terms of: 
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1. expansion coefficients of the initial fluctuation in terms of eigenmodes on the 
CH space 

2. expansion coefficients of eigenmodes on the CH space that are extended onto 
the whole pseudosphere in terms of eigenmodes on the pseudosphere. 
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A Boundary integral equation 



Here, we derive the boundary integral equation (|TT|) in section 1. For simplicity, 
we prove the formula in 3-spaces. 

First, we start with Eq.flTOD with dimesnsion M — 3. Although the integrand in 
Eq.fllOD is divergent at x = y G dQ, the integration can be regularized as follows. 
Let us draw a sphere with center y 6 dQ with small radius e and let T e be the 
outer spherical boundary and a and (3 be the internal solid angle and external solid 
angle as shown in figure 0, 

u(y)+ [ G E (x,y)V lU ^dS l - / (ViG B (x,y))« y/gdS* = 0. (63) 
The singular terms in Eq. (|63|) can be separated from non-singular terms as 



lim/ GEfcy^iUy/gdS* = / G E (x, y)V iU ^dS l 

</<9r+r e Jar 

+ lim / G E {x,y)ViU^gdS\ 

lim/ {V t G E {^y))u^dS l = / (V^x, y))u ^gdS* 

Jdr+r e J or 

+ lim / (ViGEfryfiuy/gdS* 



(64) 



If e is sufficiently small, the region enclosed by T e can be approximated as an 
Euclidean subspace. In this region, the asymptotic form of the free Green's function 
G E takes the form 

,• ^ , \ expOikd) 1 ik , nl ,, 

hmG E (x. y) = -^i = — + (65) 

where d is the Euclidean distance between x and y. Taking the spherical coordi- 
nates (e, 9, 4>) with center y, the singular terms in Eq.(|6~4|) are estimated as 



lim/ G E (x,y)VA*)V9dS l = lim - / -L^le 2 dn = 0, 
lim/ (V i G s (x,y)) M (x) v ^^ = lim / - 1 2 ^ ^(x)e 2 ^ 



^u(x), (66) 
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Figure 7: Boundary integral 

where dfl denotes the infinitesimal solid angle element. Taking the limit e — > in 
Eq.(^3]), we have the boundary integral equation in 3-spaces, 

r«(y)«(y)+/ G E (^, y )v l u^ds t - [ (v t G E (^ y ))u^ds l = o, (67) 

47r Jan Jan 

where a(y) denotes the internal solid angle at y. If the boundary is smooth at y, 
a(y) is equal to 2ir which gives the coefficients 1/2 in Eq. (p|) . Similarly, one can 
prove the formula for M = 2 and M > 3. 
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B Three-dimensional hyperbolic spaces 



The discrete subgroup T of PSL(2, C) which is the orientation-preserving isom- 
etry group of the simply-connected hyperbolic 3-space H 3 is called the Kleinian 
group. Any CH space (either manifold or orbifold) can be described as compact 
quotients M. = EI 3 /r. The classification of the Kleinian group has not been com- 
pleted. However, several procedures for constructing compact hyperbolic spaces 
are known. For further detail, see EE, |21 



The standard pseudospherical coordinates (x, 9, <p) for H 3 with curvature radius 
R are given by 

X = .Rcoshx, Xi = R sinh x sin 9 cos <fi 

X2 = R sinh x sin 9 sin 0, X3 = R sinh x cos 9 (68) 

with 

0<x<oo, 0<9 < 7T, O<0<2tt. (69) 

In these coordinates, the line element takes the form 

3 

ds 2 = Y^( dx i? = r2 W 2 + sinhx 2 (d# 2 + sin 2 ^0 2 )]. (70) 

i=0 

The Poincare representation is obtained by the transformation 

x\ = i?tanh — sin 9 cos0, X2 = -Rtanh — sin 9 sin <ft, X3 = i?tanh — cos 9, (71) 

which maps H 3 into the open ball {(xi, X2, X3) G E 3 | x\ + x\ + x\ < R 2 } called the 
Poincare ball. The line element in these coordinates takes the form 

ds2 = 4{dxl + dx\ + dxj) ^ 



R 2+ 3 



and the geodesic distance d between x and x' is given by 

2lx - x'l 2 

cosh[i?- 1 d(x, x')] = 1 + 73 ' 2W1 ' ,-., 2V (73) 



2 


x — x' : 


i 


(l"|x| 2 )(l- 


|x' 
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where | | denotes the Euclidean norm and x = i? _1 x, x' = i? _1 x'. Note that 
geodesies in the Poincare ball are either diameters or semi-circles which orthogo- 
nally intersect with the boundary of the Poincare ball. 

In Poincare coordinates, the metric is conformally flat so that the computation 
of the boundary integral equation becomes simpler. 

Another commonly used set of coordinates is obtained from the upper-half 
space representation which is defined by the transformation 

sinh x sin 9 cos <fi sinh x sin 9 sin 1 

yi = -jj ,m = d >2/3 = -p, 

D = coshx — sinh x cos 9, (74) 



which maps H 3 into the upper-half space = {(yx, t/2, 2/3) £ E 3 | y 3 > 0}. In these 
coordinates, the line element takes the form 

dg2 = R 2 (dyl + dyl + dyl) 

vi 

The geodesic distance is given by 

cosh[i?- 1 rf(y , y')] = 1 + Iy 2 ~jj - (76) 

In the upper-half space model, geodesies are either straight vertical lines or semi- 
circles orthogonal to the boundary of the upper-half space. In this coordinates 
the metric is conformally flat as in Poincare coordinates. If we represent a point 
p on the upper-half space, as a quaternion whose fourth component equals zero 
27 1, then the actions of PSL(2, C) on H 3 U C U {00} can be described by simple 



formulas, 

ap b 

^ :p -^p' = ad-bc=l, p = z + y 3 j, z = y 1 +y 2 \ ) (77) 

cp + d 

where a, b, c and d are complex numbers and 1, i and j are represented by matrices 

as, 

As p = z — y^}, the action 7 is explicitly written as 
7:i 3 UCU{oo} -> H 3 UCU{oo}, 



, . (az + b)(cz + d) + acyl y 3 \ 

^.(z(y hy2 ), y3 ) - ( \^TWTWyJ > WTWTW?J - (79) 
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If we restrict the action 7 on C U {00}, or equivalently, = 0, the action is 
described as 

. az + b 

7 : z — > z = -, ad — be = 1, z = y\ + i«2, (80) 

cz + a 

7 is called the Mobius transformation, and 7 is called the extended Mobius trans- 
formation. 

In the Klein (projective) model, the geodesies and planes are mapped into 
their Euclidean counterparts. Since the fundamental domain is enclosed by Eu- 
clidean planes in the Klein coordinates, the task of generating meshes is much 
easier than other coordinates. The transformation 

Z\ = .Rtanhxsin#cos0, Z2 = -Rtanhxsin$sin0, z% = i?tanhxcos# (81) 

can be understood as the projection of the hyperboloid (X , X%, X 2 , X 3 ) onto the 
interior of the sphere (R, z\, Z2, -23) along lines originating from the origin (0,0,0,0). 
The geodesic distance can be represented as 

cosh[iTMz, z')] = /n , =■ (82) 

where • denotes the Euclidean inner product and z = i?^ x z, z' = i? _1 z'. 

The possible values for the volume of the CH manifolds are bounded below 
and no upper bound exists. The minimal value has not yet been known, although 
Gabai et al have proved that V^o/ m j n > 0. 16668. ..-R 3 |29|. Thurston proposed a 
manifold Q2 as a candidate for the three-dimensional hyperbolic manifold of the 
minimum volume Vol(Q 2 ) = 0.98139.R 3 ||. However, Weeks f31| and indepen- 



dently, Matveev and Fomenko |32] discovered a CH manifold Qi with the smallest 
value Vol(Q\) = 0.94272i? 3 in the known CH manifolds and it is conjectured to be 
the one with the minimum volume. A computer program "SnapPea" by Jeff Weeks 
f33| has made it possible to catalog and study a large number of CH and non-CH 
spaces which include Qi,Q 2 and thousands of cusped and non-cusped hyperbolic 
3-manifolds. 
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Figure 8: Fundamental domain of the Thurston manifold which is viewed from a 
point (0,0,-1000) and a point (0,0,1000) in the Klein coordinates, respectively. 

Let us see how CH manifolds are characterized in the SnapPea. Any element of 
the discrete isometry group T which is equivalent to the fundamental group 7Ti(A4) 
can be described as a word which is a product of generators {g±, . . . , g s }, 

9 = 9mi--- 9 n m P U, nj E Z, m j = l, ...,a), g G T. (83) 

The above expression is not unique, since they are subject to a set of relations, 
each of which takes the form, 

J (j,kjeZ, l, = l, ...,s), (84) 

where I denotes the identity. Note that different expression of g is possible by 
choosing different generators. In the case of Thurston's manifold Q2, T has a 
simple presentation, 

r = {a, b : ab 3 aba~ 2 b, ab~ 1 a~ 1 baba~ 1 b~ 1 ab~ 1 }, (85) 

where a and b are generators and words in the parenthesis are equal to identities. 

This representation is simple for describing F but not convenient for describing the 

fundamental domain. Choosing a coordinate system centered at a point of locally 

maximum of the injectivity radiusf], generators which define the face identification 

7 The injectivity radius of a point p is equal to half the length of the shortest periodic geodesic 
on p. 
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maps in the pseudospherical coordinates can be described by 8 matrices (see ap- 
pendix C) , which implies that the number of the faces on the boundary of the fun- 
damental domain is sixteen. For instance, the center (X , Xi, X 2 , X 3 ) = (1, 0, 0, 0) 
is moved to a point (1.63032, -0.5657, 0.993147, -0.592943) by T\. The fundamen- 
tal domain of the Thurston manifold can be computed from these 8 matrices. First, 
let us make 16 copies of the basepoint (1, 0, 0, 0) that are obtained by multiplying 
the basepoint by 8 matrices and their inverse matrices. Next, the basepoint and 
the 16 copies are connected with 16 geodesic segments. If one put planes on the 
equidistant points on the segments at right angle, then one obtains the fundamen- 
tal domain enclosed by the 16 planes. The face identifications are shown in figure 
H in which each color of the faces corresponds to one of the identification maps. 
For instance, a point on the face with red color in the left figure is identified with 
the corresponding point on the face with the same color in the right figure by T\. 
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C Table of matrices 



In the Minkowski coordinates (t,x,y,z), the 8 generators which define the funda- 
mental domain of the Thurston manifold with the basepoint (1, 0, 0, 0) are described 
by the following 8 matrices, 



7i 



To 



T, = 



/ 1.630319018827 
-0.565700099811 
0.993147485366 
\ -0.592943144435 
/ 1.630319018827 
-0.725744142114 
-0.749088496222 
\ 0.755050970573 
/ 1.630319018827 
-0.493618131693 
0.044016981085 
\ -1.188420695119 
/ 1.630319018827 
0.990384655015 
-0.430529153566 
\ -0.701229624437 
/ 1.630319018827 
-1.169223329346 
-0.48230025187 
\ -0.241336645178 
/ 1.630319018827 
-0.334091294985 
-1.230309911716 
\ -0.180722524595 



-0.823154893134 
0.926953540114 
-0.896311100897 



0.949658853916 
0.03235859132 
1.03061185404 



-0.122341828367 -0.915775117974 
-0.404520784012 -0.353128575811 



1.032680902406 
0.153460453167 
0.271398430022 

0.687449077 
0.144914737528 
0.704808911001 
-0.977154210597 



0.430856766432 
-0.490988659119 



-0.280185985083 \ 
0.678031458683 
0.34757379707 
0.705669882356 / 

-1.1 70300494487 \ 
0.524058355949 
1.138645511607 



-0.835459379318 -0.893561685534/ 
-0.510375368061 -0.961702060597\ 



1.075867816726 
0.053046648969 
0.874394274525 / 



-0.255278318146 
0.708779194553 
0.832435012465 
-0.076636501663 -1.186821280266 -0.493479683947 \ 
0.188025082259 -1.393877403859 0.051128429897 
0.444575008386 
0.517690637565 



0.973085065179 
-0.153704839725 



-0.422808564099 -0.669242281726 
0.245239342961 1.092969619938 



1.00284583704 
-0.336043112815 
-1.258320677027 
0.339168543325 
1.513335320951 
0.422080361796 



0.202024429006 
1.09548811596 / 
1.015523406733 \ 
-1.054683966146 
-0.146468308544 -0.453277833869 



-0.481507711397 
-0.085154694975 

0.989853092692 
-0.161421439076 

0.03721956436 



-0.84467077995 / 
-0.259456845811 \ 
0.129508870322 
0.444321735047 
-0.92364684978 / 
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/ 1.815210205228 
-1.394305392881 
0.587810577053 

V 0.073343616291 
/ 1.815210205228 

0.093635040681 
-0.930634634266 

V —1.191696163383 



-1.011761093044 
1.36137877709 
0.152666848297 

-0.383407313596 

-0.288909259155 
0.602143185543 
0.833203272177 

-0.163292533092 



-1.117355434627 
0.985732513048 
-1.125213789716 
0.103481918443 
0.774531377414 
-0.536756936203 
-0.104192832263 
-1.1405852444 



0.15114368155 \ 
-0.345060526033 
-0.236871347696 
-0.920689747139 / 
-1.269496228642 \ 
-0.598400448183 
1.077495699961 
1.045246666925 / 
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